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Abstract. The improvement of simulations of QCD with dynamical Wil- 
son fermions by combining the Hybrid Monte Carlo algorithm with paral- 
lel tempering is studied. As an indicator for decorrelation the topological 
charge is used. 



1. Introduction 

Decorrelation of the topological charge in Hybrid Monte Carlo (HMC) sim- 
ulations of QCD with dynamical fermions is a long standing problem. For 
staggered fermions an insufficient tunneling rate of the topological charge 
Qt has been observed [1, 2]. For Wilson fermions the tunneling rate is ad- 
equate in many cases [3, 4]. However on large lattices and for large values 
of k near the chiral limit the distribution of Qt is not symmetric even after 
more than 3000 trajectories (see Figure 1 of [3] and similar observations by 
CP-PACS [5]). 

It has also been observed that sensitive observables like the r( corre- 
lator are Qt dependent [10]. Thus it appears to be important to look for 
simulation methods that give good distributions of Qt- 

The idea of parallel tempering is to improve transitions in parameter 
regions where tunneling is suppressed by opening ways through parameter 
regions with little suppression. In QCD the method has been applied suc- 
cessfully for staggered fermions [6] . In [7] parallel tempering has been used 
to simulate QCD with 0(a)-improved Wilson fermions without finding any 
gain, however, with only two ensembles which does not take advantage of 
the main idea of the method. 

Here parallel tempering is used in conjunction with HMC to simulate 
QCD with (standard) Wilson fermions. The gain achieved is demonstrated 
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by studying time series and histograms of the topological charge and by 
comparing statistical errors of the topological susceptibility (Qt)- 



2. Parallel Tempering 

In standard Monte Carlo simulations one deals with one parameter set A 
and generates a sequence of configurations C. The set A here includes (3, 
k, the leapfrog time step and the number of time steps. C comprises the 
gauge field and the pseudo fermion field. 

In the parallel tempering approach [8, 9] one simulates N ensembles 
(Aj;Cj), i = 1, . . . , N in a single combined run. Two steps alternate: (a) 
update of N configurations in the standard way, (b) exchange of configu- 
rations by swapping pairs. Swapping of a pair of configurations means 

(l\ r wi r w _^ / (( A ^ C j)' ( x j> c i))> ifaccepted m 

(,(A, Ci), (Aj, Cj)) -> | i{x .. c . )AXj . Cj) ^ else (i) 

with the Metropolis acceptance condition 

P swap (i,j) = mm (l, e~ AH ) , (2) 

AH = H\.(Ci) + H x .(Cj) - H Xi (Cj) - H Xj (C t ). (3) 

Since after swapping both ensembles remain in equilibrium, the swapping 
sequence can be freely chosen. In order to achieve a high swap acceptance 
rate one will only try to swap (/?, «;)-pairs that are close together. If the 
chosen (/?, «;)-values lie on a curve in the (/?, «)-plane there are three obvious 
choices for the swapping sequence of neighboring K)-pairs. One can step 
through the curve in either direction or swap randomly. It has turned out 
that it is advantageous to step along such a curve in the direction from 
high to low tunneling rates of Qt- 



3. Simulation Details 

The standard Wilson action for the gauge and the fermion fields was used. 
The lattice size was 8 4 . The HMC program applied the standard conjugate 
gradient inverter with even/odd preconditioning. The trajectory length was 
always 1. The time steps were adjusted to get acceptance rates of about 
70%. In all cases 1000 trajectories were generated (plus 50-100 trajectories 
for thermalization) . 

Qt was measured by the field-theoretic method after 50 cooling steps 
of Cabibbo-Marinari type. This method gives close to integer values which 
were rounded to the nearest integers. (Note that the results presented in 
[11] were obtained without rounding.) 
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Statistical errors were obtained by binning, i.e., the values given are the 
maximal errors calculated after blocking the data into bins of sizes 10, 20, 50 
and 100. 

4. Results 

Several tempered HMC simulations were run in the quenched approxima- 
tion (tempering in (3) and with dynamical fermions (tempering in k, at fixed 
(3 = 5.5 and (3 = 5.6). For comparison also standard HMC simulations have 
been performed. 

Figures 1 and 2 show typical comparisons of time series and histograms 
of Qt- One sees that with tempering considerably more topologically non- 
trivial configurations occur and that the histograms of Qt become in general 
more symmetrical and broader. 

In standard runs Qt frequently stayed for quite some time near 1 or 
near —1, while with tempering this never occurred. The standard run at 
k = 0.156 shown in Figure 2, where Qt gets trapped in this way for about 
200 trajectories, provides an example of this. Such observations have also 
been made on large lattices [3, 5]. 

While a correlation analysis cannot be carried out with the given size 
of samples, some quantitative account of the improvement by tempering is 
possible using the mean of the absolute change of Qt, called mobility in [3], 

D i = 7j— E \Qt(*)-Qt(i-i)\ • (4) 

Results for D\ are given in Tables 1 and 2. If \Qt{i) — Qt(i — 1)\ < 1 for all 
trajectories then l/D\ is the HMC time between topological events. Since 
that condition holds in most of the cases presented here one gets an idea 
of the quantitative improvement by tempering. 

Another quantitative estimate of improvement comes from the statisti- 
cal errors of (Qt)- The fact that statistical errors decrease with the square 
of HMC time provides a second quantitative criterion for the speed-up of 
a simulation. 

Quantitative results at (3 = 5.5 are summarized in Table 1. From the 
ratios of mobilities and squared ratios of errors of susceptibilities one ob- 
tains speed-ups between 2 (ratio of D\ at k = 0.158) and 16 (squared ratio 
of the errors of (Ql) at n = 0.160). This is a considerable gain, especially 
if runs at several values of k need to be done, what is usually the case. 

At (3 = 5.6 tempering looks even better in the sense that the standard 
HMC runs do not really resolve the topological properties for k > 0.156 
(see Table 2 and Figures 2 and 3). 
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Figure 1. Comparison of time series and histograms of Qt obtained from standard and 
tempered HMC on the 8 4 lattice at (3 = 5.5. In the tempered run 5 ensembles were used, 
0.158 < k < 0.160 and Ak = 0.0005. The swap acceptance rate was about 56%. 



TABLE 1. Mobilities D\ and topological susceptibilities 
(Qt) f° r the plots shown in Figure 1. 





Standard HMC 


Tempered HMC 


K 


D 1 


<Q?> 


Di (Qt) 


0.158 


0.171(35) 


0.51(19) 


0.398(53) 0.49(8) 


0.159 


0.058(20) 


0.12(5) 


0.248(40) 0.20(5) 


0.160 


0.012(8) 


0.030(27) 


0.056(13) 0.031(7) 



In the following the choice of K-values at j3 = 5.6 is motivated. The 
run with 21 ensembles can be considered as a reference run. In a large 
scale simulation one would want to use less ensembles. The run with 6 
ensembles demonstrates that comparable speed-up can be achieved with 
a smaller number of ensembles. The run with 7 ensembles covers exactly 
the parameter range investigated by SES AM [3] . It was mainly done to get 
estimates for the swap acceptance rate on larger lattices for Ak = 0.00025 
(see section 5). 

It is interesting to compare the runs with 6 and 7 ensembles. In the run 
with 6 ensembles the mobility is higher. This reflects the main idea of the 
tempering method which is to connect areas of low tunneling rates with 
areas of high tunneling rates. 




Figure 2. Comparison of time series and histograms of Qt obtained from standard and 
tempered HMC on the 8 4 lattice at (i = 5.6. The corresponding quantitative results can 
be found in Table 2. 
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TABLE 2. Mobilities D\ and topological susceptibilities (Qt) on the 8 4 lattice at f) = 5.6. 
The swap acceptance rates achieved were about 82% for Ak — 0.00025 and about 63% for 
Ak = 0.0005. 





Standard HMC 




Tempered HMC 








7 ensembles 


6 ensembles 


21 ensembles 






0.156 < k < 0.1575 0.155 < k < 0.1575 


0.15 < k < 0.16 






Ak = 0.00025 


Ak = 0.0005 


Zlft = 0.0000 


K 










0.1500 


0.325(39) 






0.764(46) 


n i ^9n 


fl 1 7A( r >Q\ 






n 7 <i i t i(' c in , i 

u. t ooyou j 


A 1 C KO 

O.looO 


0.031(11 ) 




/~\ -t Sit—?/ A C\\ 

0.167(42) 


U.L6Z(6Z ) 


fl 1 






0.176(43) 


fl 1 1 ^n"i 

U. 1 lO^OU J 


U. 10DU 


U.UoU^l 1 ) 


0.064(27) 


0.108(36) 




n i 




0.040(15) 


0.102(28) 


f) 074^99"! 

U.U / ^iyLL J 


U. 10 / u 


U.UUZ^Z J 


U.Uzz^o J 


n n^fi/i o\ 
u.uoo^iy ) 




n i 




0.004(3) 


0.044(14) 












U.Ulu^O J 


0.1600 









n 



K 






(Qt) 




0.1500 


0.77(14) 






0.993(92) 


0.1520 


0.58(13) 






0.707(57) 


0.1550 


0.056(29) 




0.144(40) 


0.085(21) 


0.1555 






0.100(25) 


0.071(18) 


0.1560 


0.134(83) 


0.044(20) 


0.062(23) 


0.052(14) 


0.1565 




0.020(8) 


0.055(17) 


0.040(14) 


0.1570 


0.004(4) 


0.011(4) 


0.037(11) 


0.028(9) 


0.1575 




0.002(1) 


0.030(11) 


0.017(6) 


0.1580 


0.004(4) 






0.008(4) 


0.1600 













5. Going to larger Lattices 

With regard to large scale simulations of QCD performance predictions are 
needed. One potential problem of the tempering method has been stressed 
in [7] , namely the decrease of the swap acceptance rate (A) with the lattice 
volume. In [7] it has been checked that the relation [12] 

<A) = erfcQ v /(W) (5) 

is valid for a large range of (AH). Relation (5) also holds in all simulation 
done in this work. 
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Figure 3. Comparison of topological susceptibilities on the 8 4 lattice at f3 = 5.6. The 
plot shows results from standard HMC (N = 1) and the tempered HMC run with N = 21 
ensembles. 

Because {AH) scales linearly with the lattice volume V, relation (5) 
allows one to predict (A) by inserting values measured on the 8 4 lattice. 
Table 3 lists predictions using values of (A) from the runs shown in Table 2. 
Some caution is necessary with these predictions because on the 8 4 lattice 
at (3 = 5.6 and 0.15 < k < 0.16 the finite temperature phase transition [13] 
is crossed. 

TABLE 3. 

Ak V (A) 

0.0005 8 4 63% 
16 3 x 32 0.6% 

0.00025 8 4 82% 
16 3 x 32 20% 
24 3 x 48 0.4% 



Indeed more and more ensembles will be needed on larger lattices if one 
wants to keep (A) and the parameter range constant. However it is an open 
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question which effect is stronger, the decrease of {A) or the slowing down 
of tunneling between topological sectors. The hope is that the need to take 
more ensembles more than compensates the slowing down of tunneling. 

6. Conclusions 

On the 8 lattice parallel tempering considerably enhances tunneling be- 
tween different sectors of topological charge and generates samples with 
more symmetrical charge distributions than can be obtained by standard 
HMC. The histograms also get slightly broader or even become nontrivial 
thanks to this technique. 

The enhancement of tunneling indicates an improvement of decorrela- 
tion also for other observables. More satisfactory histograms are important 
for topologically sensitive quantities. Both of these features make paral- 
lel tempering an attractive method for large-scale QCD simulations. The 
method is particularly economical when several parameter values have to 
be studied anyway. 

A potential problem is that for a given parameter set the swap ac- 
ceptance rate (2) decreases for increasing lattice volume [7]. To settle the 
question whether on larger lattices the need for increasing the number of en- 
sembles is compensated by improved tunneling between topological sectors 
this study will be continued on larger lattices. 
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